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The vortex-vortex interaction potential in bulk superconductors is calculated within the Ginzburg- 
Landau (GL) theory and is obtained from a numerical solution of a set of two coupled non-linear 
GL differential equations for the vector potential and the superconducting order parameter, where 
the merger of vortices into a giant vortex is allowed. Further, the interaction potentials between a 
vortex and a giant vortex and between a vortex and an antivortex are obtained for both type-I and 
type-II superconductors. Our numerical results agree asymptotically with the analytical expressions 
for large inter-vortex separations which are available in the literature. We propose new empirical 
expressions valid over the full interaction range, which are fitted to our numerical data for different 
values of the GL parameter. 
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I. INTRODUCTION 



The interaction potential between vortices has been 
an important study subject for many years. In 1971, 
Kramer— used the asymptotic behavior of the vortex 
fields for large distances in the Abelian Higgs model to 
obtain an analytical expression for the vortex- vortex in- 
teraction potential, which is given by modified Bessel 
functions. This potential is attractive (repulsive) for 
type-I (type-II) systems, i.e., for a superconductor 
Ginzburg-Landau parameter with k = A/£ < l/v2 
(> 1/V2)) where A is the penetration depth and £ is 
the coherence length. Moreover, this expression leads to 
a constant interaction energy as a function of the separa- 
tion between vortices for the critical value n = l/\/2 (also 
called the Bogomol'nyi point), implying that vortices do 
not interact in this regime. A detailed analysis of the 
vortex-vortex (V-V) interaction was given later by Ja- 
cobs and Rebbi, 2 who constructed a variational function 
describing two separate vortices and obtained the vari- 
ational parameters by minimizing the free energy. This 
variational function was able to model: i) the deforma- 
tion of the vortex core when the vortices are brought close 
to each other, and ii) the formation of a giant vortex^r— 
when vortices are superimposed on each other. 

Thereafter, several works studied different aspects of 
the interaction between vortices in superconductors. For 
example, Brandt used the asymptotic expression for the 
interaction potential to study the elastic properties of 
flux-line lattices in type-II superconductors. Speight^ 
derived the V-V interaction from a linear field theory 
described by a Lagrangian of two singular point sources 
placed at the vortex centers. MacKenzie et alM- used 
the linear field theory proposed by Speight to obtain 
the interaction between separated strings in a model 
with two order parameters, which may be relevant for 
superconducting cosmic strings^ the SO (5) model of 



high-temperature superconductivity and solitons in non- 
linear optics. Eventually, all these analytical approxi- 
mations lead to expressions which are the same, or at 
least very similar, to the one derived by Kramer for 
large separation between vortices. Other models were 
presented by Mohamed et al. 10 , who used a perturba- 
tive approach to calculate the V-V interaction for super- 
conductors with k ~ l/v2 when the Ginzburg-Landau 
(GL) theory is extended to low temperatures, and by 
Hernandez and Lopezjii who used a variational approach 
based on the Clem trial function^ to calculate the force 
between vortices. Auzzi et al~2- showed recently that 
for non- Abelian vortex interactions, there are two other 
regimes besides the well-known type-I and type-II: de- 
pending on the relative orientation, the interaction po- 
tential can show attractive and repulsive regions for the 
same system. Similar behavior can be achieved in a two- 
band superconductor! 14 i 15 

In the present work, we solve numerically the set of 
Ginzburg-Landau (GL) differential equations for two vor- 
tices fixed at a certain distance from each other. We 
are able to recover the V-V potential obtained by Ja- 
cobs and Rebbi 2 . and the asymptotic behaviors predicted 
by Kramer—. The results are generalized to the case of 
vortex-giant vortex (V-GV) and vortex-antivortex (V- 
AV) interactions. We parameterize all obtained forces 
for future use in e. g. molecular dynamics simulations of 
the vortex matter. 

The remainder of this paper is organized as follows: in 
Sec. II, we review the asymptotic behavior of the mag- 
netic field and order parameter of the single static vortex 
in the GL theory and expand the analytical expression 
suggested by Kramer to study the interaction between a 
single vortex and another vortex, an anti vortex or a gi- 
ant vortex, in the limit of large separation between them. 
In Sec. Ill we discuss the set of coupled non-linear differ- 
ential equations that describe these interactions, which 
are valid for arbitrary separation between vortices, and 



thus, also for small distances where deformations of the 
interacting vortex cores are important. This set of equa- 
tions is solved numerically without any approximations 
for arbitrary values of k, and the results for the V-V, 
V-GV and V-AV interactions are shown in Sec. IV and 
compared to the analytical expressions obtained in Sec. 
II. For each of these cases, a fitting function for the inter- 
action is proposed and the fitting parameters are given. 
Our results are summarized in Sec. V. 



II. ASYMPTOTICS OF THE INTER- VORTEX 
POTENTIAL 

Let us start with the expression for the free energy 
in the GL theory, or equivalently, the potential in the 
Abelian Higgs model 2 -^: 



where the functional T is given by 
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(2) 
In this expression, tp is the order parameter (a complex 

scalar held), A is the electromagnetic vector (gauge) po- 
tential and a and j3 are phenomenological parameters, 
which are related to the two characteristic lengths in a 
superconductor: the coherence length £ = %/ ' \fAma and 
the penetration depth A = c/e(^ymf3/8^ra). It is conve- 
nient to define /i = v2k, so that when /i < 1 (/i > 1) 
we are in the type-I (type-II) regime (« = l/\/2 leads to 
(j, = I). 2 - If we define A = (2e/hc)A and * = ^/jT/^^/X, 
the expression for the energy functional can be rewritten 
in dimensionless units as 
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where the distances are scaled by the magnetic field pen- 
etration depth A, energy by Eq — /?/2a 2 £ 2 and the force 
by fl Q = (3/2a 2 £ 2 \. 

The lowest energy configuration of the system is found 
by minimizing E with respect to \& and the vector po- 
tential A. The standard way of minimizing a functional 
is by the Euler-Lagrange equations, which in the case of 
Eq. fll]) are 
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Equations (14a)) and (|4bl) result in 
Ginzburg-Landau equations. 16 
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the well-known 



To find the V-V interaction potential, one has to con- 
trol the localization and the winding number (also called 
vorticity) of the vortices. As mentioned earlier, previ- 
ous theoretical works propose a way to fix the vortices 
and obtain the interaction potential, based on four steps: 
i) fixing a circular phase change of 2tt around each vor- 
tex, ii) finding the asymptotic behaviors of the vector po- 
tential and the amplitude far from each of the vortices^ 
or, equivalently, solving numerically the differential equa- 
tions for these variables? 2 - iii) constructing a superposi- 
tion ansatz for ^ and A which describes the double vortex 
structure, and iv) using the latter ansatz in Eqs. ^ and 
(|3]) to obtain the energy for a given separation between 
vortices. 

If we have the asymptotics of the vortices, obtained in 
step ii), the integral in Eq. (fTJ) can be solved analyti- 
cally, giving an analytical expression for the energy as a 
function of the separation between vortices. In order to 
find the analytical expression for the asymptotic interac- 
tion potential between a vortex and another vortex, an 
antivortex or a giant- vortex, we start from the sequence 
listed above: a circular phase is fixed around each vortex 
by assuming ^(r, 9) = f(r)e m6 , where (r, 0) are polar 
coordinates with the origin in the center of the vortex, 
n is its winding number and f(r) is the amplitude of its 
order parameter, which is assumed to be circularly sym- 
metric around the vortex center. Considering the gauge 
A = na(r)0/r, the Euler-Lagrange equations (0a| and 
(|4b[) for a single vortex read 
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with /(oo) = a(oo) = 1. As is well known 2 -, the winding 
number n also determines the number of zeros of the 
vortex field f(r) and, because of the circular symmetry, 
these zeros must be degenerate at r = 0. Substituting the 
auxiliary functions a{r) = f(r) — 1 and Q(r) = a(r) — 1 
in Eqs. (|5ap and (|5b[) . we can suppress the high order 
terms in the remaining differential equations when r — ¥ 
co, as er(oo) = Q(oo) = 0, which leads to the following 
equations, valid in the asymptotic limit: 
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Equations (I6al) and (|6bl) are easily identified as modified 
Bessel equations and their solutions are a(r) = 7 1 _ftr (/ ir ) 
and Q(r) = "f2fKi(r) where 71 and 72 are coefficients 
to be determined. For example, after solving Eqs. (|5a[) 



and (I5bl) numerically, one can obtain these coefficients by 
fitting /(r) = l+jiK (jtr) and A = n{l + ^ 2 rKi(r))0 /r 
to the results obtained by the numerical procedure. 

Several different procedures can be followed to extract 
the analytical expression for the interaction potential 
from these asymptotic functions. As an example, Betten- 
court and Rivers 9 suggested that one can substitute the 
superposition ansatz ^l{r,r\,rz) = \&(|r — ri|)\fr(|r — r 2 |) 
and Ae(r,ri,r 2 ) = A$(\r — ri\) + A$(\r — r2\), for vortices 
centered at f[ and rjj, in the energy functional in Eq. (TTJ) 
and keep only the terms that are linear in the fields for 
each vortex. After some calculations described in more 
detail in Refi^, one obtains 

E int (d) = 2 [mrrf ] l ( 2 ] K a {d) - 7 j 1) 7 f ) ^o(Mrf)] , (7) 

where d is the separation between vortices, Hi is the wind- 
ing number and 7} and 72 are the fitting coefficients 
for the vortex i in position f\. We point out that in the 
paper by Bettencourt and Rivers, the expression for Ei nt 
is slightly different from Eq. ([7]) because they assumed 

that 7} = |?ij|72 1 which is valid only in the critical 
coupling regime for fi = 1, as stressed by Speight 7 - and 
verified by BogomoPnyi equations for this regime. The 
same expression for Ei nt (d) was found by Kramer by a 
perturbational approach 1 and can also be obtained by 
considering point sources in a linearized field theory 7 ' 8 . 

It can be easily seen that the potential in Eq. (J7J) is 
consistent with the fact that for fi < 1 (> 1), the V-V 
interaction is attractive (repulsive) for large separation, 
leading to a type-I (type-II) superconducting behavior. 
This statement is also valid for vortex-giant vortex in- 
teractions, as Eq. ([7]) still holds for this case, where the 
interaction potential is obtained just by setting m = 1 

and Tin > 1 and finding the fitting coefficients 7} and 7 2 
for this case. However, for a V-AV interaction, n\ni is al- 
ways negative, leading to an attractive potential Ei nt (d) 
for any value of k. Furthermore, Eq. (J7J shows that 
the interaction between an antivortex and a giant vor- 
tex is always attractive as well. This can be understood 
through a heuristic argument: when a vortex and an an- 
tivortex are far from each other, the energy of the system 
is non-zero, as it is the sum of the energies of one vor- 
tex and one antivortex; on the other hand, when they 
approach each other they should annihilate, giving zero 
energy. Hence, at least at some distance, the energy of 
the V-AV pair must decrease as d approaches zero and, as 
a result, the interaction potential is attractive. This re- 
sult is in contradiction with Ref^, where it was claimed 
that in type-I superconductors the interaction of a V- 
AV pair is repulsive, and which was used to explain the 
existence of a stable V-AV molecule in mesoscopic su- 
perconducting triangles 18,19 . As follows from our theory, 
the V-AV interaction is always attractive and should lead 
to a disfavored V-AV molecule in type-I superconducting 
polygons, because in that case the V-V interaction is also 



One more question arises from the conjecture of Ref. 17 : 
if vortices attract (repel) each other in type-I (type-II) 
superconductors, whereas exactly the opposite occurs for 
V-AV interactions, what one would expect in a critically 
coupled system? In this regime, i.e., when fj, = 1, vortices 
do not interact; should we expect the same for V-AV? 
The answer to this question is provided by the Bogo- 
mol'nyi equations 22 
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For the single vortex ansatz chosen before, ^(r, 9 
f(r)e m6 and A = na(r)6/r, these equations read 2 ^ 
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Substituting the latter formulae in Eq. (j9]) and ne- 
glecting higher-order terms in a(r) and Q(r) yields 
nQ = —sgn(n)rjr; or, using the asymptotic forms of 
these functions, WYzrKiir) = —sgn(n)"firdKo(r)/dr — 
sgn(n)jirKi(r) => 71 = sgn(n)nrf2- Substituting this 
expression for 7} in Eq. ([Jj shows that the interac- 
tion potential for V-AV in the critical coupling regime 
is Ei n t(d) = 4rti72 "-272 Ko(d) : which is still attrac- 
tive since n\ni < 0. Hence, unlike vortex-vortex pairs, a 
V-AV pair exhibits an attractive interaction even in the 
critical case of fi = 1. 

It should be mentioned that, as 71 7^ 72 for [i ^ 1, 
the interaction potential given by Eq. (J7J) may diverge 
at small distances for some values of \i. This signals the 
breakdown of this analytical expression for small d. Actu- 
ally, for small separation d, the superposition ansatz pro- 
posed by Bettencourt and Rivers and used in the present 
section also fails, since it does not take into account ei- 
ther the spatial deformation of the fields, or the possi- 
bility of the formation of giant vortices. 2 Hence, the an- 
alytical expression for the interaction potential between 
vortices has a validity restricted to large d. Nevertheless, 
when the interaction potential is calculated by numeri- 
cal means, taking into account all the features mentioned 
above, the result shows very good agreement with Eq. (J7J 
for d larger than a critical separation d c which depends 
on (a, as will be shown in the next section. 

The interaction force fi can be obtained by taking the 
derivative of the energies with respect to the distance 
between the two vortices. The results for /i = 1 in the V- 
V and V-GV cases clearly give 0(d) = 0. The analytical 
expression for the force for large vortex- vortex separation 
can be easily derived from Eq. (J7J as 

n(d) = 2 In^^^K^d) - 7 1 1) 7i (2) M #iM1 • (10) 



III. GL EQUATIONS FOR FIXED VORTICES 



energy functional F for fixed position of vortices as 



We recall step i) in the procedure described in the pre- 
vious section for obtaining the V-V interaction, which is 
fixing a circular phase change around each vortex. In the 
present section, we derive the Euler-Lagrange equations 
from Eqs. (l4al) and (|4b[) with the constraints imposed by 
step i). 

In the paper by Jacobs and Rebbi 2 -, the authors fixed 
the phase for a single vortex and obtained 'modified' GL 
equations, given by Eq. (2.18) of their paper or, equiva- 
lently, Eqs. (|5a|) and (|5b|) of the present work, as well as 
Eq. (7) of the paper by Babaev and Speight^. Although 
they presented these equations in their paper, Jacobs and 
Rebbi did not solve them directly, but used instead vari- 
ational functions for f(r) and a(r) and minimized the 
energy E without solving the differential equations. 

For the two vortex system, Jacobs and Rebbi made 
a different ansatz, ^ = exp[i#i] exp[i02]f(r,9), where d\ 
and 9 2 are azimuthal angles around each vortex position, 
/ is a real function which is not necessarily circularly 
symmetric and is zero at the position of each vortex, 
and the winding numbers were chosen as 1 for each vor- 
tex. Having fixed the positions and vorticities, they just 
needed to find / and A that minimize E. As before, in- 
stead of deriving Euler-Lagrange (differential) equations, 
they used a variational procedure, considering trial func- 
tions that account for the deformation of the vortices 
towards the formation of the giant vortex. 

The results of Jacobs and Rebbi are rather accurate 
and the advantage of their approach is that many terms 
of the variational function can be integrated analytically. 
However, the variational procedure involves many param- 
eters, it is a very long analytical calculation and their 
trial function is not the most general function: if one 
wants to solve the problem for a V-AV or a V-GV pair, 
the trial function has to be modified and consequently 
also the analytical integrals in the variational procedure. 

In order to obtain the inter-vortex potential, our ap- 
proach starts from the ansatz for two vortices \E r = 
e mi 1 e m2 2 f(x, y), where we control the vorticity n\ and 
n 2 of each vortex. We further rewrite e mj6j in Cartesian 
coordinates: 
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Notice that although X and Y seem to have no physical 
meaning, they can be related to the gauge proposed by 
Jacobs and Rebbi for the vector potential, A = na(r)0/r, 
or A = (~nsm0/r,ncosd/r,O) = (—ny/r 2 ,nx/r 2 ,0) in 
Cartesian coordinates, which was shown to be compati- 
ble with the symmetry of the Euler-Lagrange equations, 
leading to Eqs. (|5aj) and (|5b)) for a single vortex. How- 
ever, as we are looking for general differential equations 
for two vortices, we will not make any a priori choice of 
gauge for the vector potential. 

In previous works,— ^ the Euler-Lagrange equations 
were not explicitly derived. We derived the Euler- 
Lagrange equations for the present problem setting, 
which are given by 
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where the unitary angular vectors around each vortex can 
be rewritten as 8j = {—yj/rj,Xj/rj,0). One can even 
verify, after some manipulations of the equations, that 
inserting n 2 = in these equations leads to Eq. (2.18) 
in the paper by Jacobs and Rebbi, which is the equation 
for a single vortex, but in Cartesian coordinates. Solving 
Eqs. (|13a|) and (|13bl) is equally demanding as solving the 
common GL equations, which was done in many works 
in the literature ^ 20 ' 21 
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where fj = (xj , y 3 - , 0) is the in-plane position vector with 
origin in the center of the vortex j. As we will study the 
case for two vortices separated by a distance d, we take 
n = (x-d/2,y,Q) and r 2 = (x + d/2,y,0). 

Next, we substitute this ansatz into Eq. ^ to get the 



IV. 



NUMERICAL RESULTS AND FITTING 
FUNCTIONS 



We solved Eqs. (I13a[) and (|13bl) numerically using 
the finite difference technique and a relaxation method 
suitable for non-linear differential equations/ 2 ^ The two- 
dimensional system is divided in a uniform square 
601x601 grid with total dimensions 60A x 60A. The sin- 
gularities in the amplitude of the order parameter appear 
naturally in the center of each vortex position, as a con- 
sequence of the fixed circular phase e mi9i defined around 



Vortex-vortex interaction 



0.12 


1 i ' i ' i ' i ' i ' i ' i 
.(a) f\ 




j \ H=l-2 


0.08 


; \ - - n=i-4 

- ; . . • — n-i.8 

1/ '\ U = 2.2 

i: \ 


,-,0.04 
G 0.00 

-0.01 




-0.02 


v/" ^ .' s - - n = o-3- 

... - . ji-0.4 




H-0.6' 


-0.03 


H = 0.8_ 






FIG. 1: (Color online) Numerically obtained vortex- vortex 
interaction force Q as a function of the separation d between 
vortices, for several values of /x — V2k in the (a) type-II and 
(b) type-I regimes, (c) Critical separation d c (open squares, 
right scale) and extremum Q ma x (open triangles, left scale), 
which correspond respectively to the position and amplitude 
of the peak in the force, as a function of the GL parameter. 
The fitting functions for d c and flmax are plotted by the solid 
curves. 



each vortex i, which guarantees the existence of zeros of 
the order parameter in the center of the vortices.— The re- 
sults for V-V, V-GV and V-AV interactions are presented 
separately in the following subsections. Analytical fitting 
functions will be proposed for the numerically obtained 
curves, where the fitting error is defined by the variance 2 ^. 
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[G(n)-G fit (n)Y 
(N - N p ) 
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where G(n) is the numerical data set, Gfit{n) is the an- 
alytical fitting function, N is the length of the data set 
and N p is the number of variational parameters of the 
fitting function^ 



The numerical results for the V-V interaction force are 
shown in Fig. [TJ for several values of \x in the type-II 
(a) and type-I (b) regimes. Notice that for /x ss vor- 
tices should not interact and the force vanishes, but for 
< (i < 1 they attract (type-I regime) and the force 
is negative in this case. However, at the critical point 
/x = 1 the force vanishes again. Hence, in the type-I 
case, two different regimes can be identified: one where 
the force increases from zero, as /x increases from zero, 
and the other when the force decreases back to zero, 
as fj, approaches 1. This can be seen from Fig. [He), 
where extremum of the force peak Q m ax (open trian- 
gles, left scale) increases with /x for small /x, but for 
0.6 < /x < 1 the peak decreases with /x, approaching 
zero when \i — 1. The numerical results can be fitted to 
OnazO) = 0.096lAi(ii - 1)/(1 + 0.2863/i) 1341 , which is 
shown by the red curve in Fig. Hfc). From Figs. QIa,b) 
we see that the force exhibits a maximum at some crit- 
ical separation d Cl which depends on the GL parameter 
/i = y/2n. The critical separation d c is also shown in Fig. 
[IJc) (open squares, right scale) as a function of /i. The 
fit of this curve, given by d c = 22.203(1 + 10.504m) -0 - 774 
(with estimated variance v sw 0.3%), is shown by the 
solid curve in Fig. [He), which suggests that the criti- 
cal separation for the V-V interaction approaches zero 
in the extreme type-II situation (/x —¥ oo). An effec- 
tive extreme type-II scenario can also be achieved in a 
superconducting film of thickness w <C A, where the ef- 
fective penetration depth is A = \ 2 /w, and for which 
analytical expressions for the V-V interaction force were 
proposed by Pearl2£, and later by Brandt^ 7 -. However, 
in the case of superconducting films, the V-V force de- 
cays monotonically as 1/d 2 , while in the present case of 
a bulk superconductor the decay is exponential. Hence, 
although both situations can be considered as extreme 
type-II limits, our results for bulk superconductors with 
fi — > oo are quantitatively different from those for thin 
superconducting films. 

The importance of solving Eqs. (I13ap and (I13b[) nu- 
merically for two separate vortices lies in the possibility 
of obtaining the interaction force between vortices even 
in the small separation limit, which cannot be described 
by the asymptotic functions given in the literature 9 and 
described by Eq. (TTUl) . However, since solving these equa- 
tions is generally not an easy task, we attempt to pro- 
pose here an analytical expression that possesses all the 
features of the numerically obtained force as a function 
of the vortex separation. Such an analytical expression 
can be helpful e.g. for numerical modelling of vortex 
structures by means of molecular dynamics (MD) sim- 
ulations, where vortices are considered as point parti- 
cles. Up to now, in those MD simulations one typically 
uses the asymptotic analytical expressions for the V-V 
interaction which are only valid in the large separation 
limitj 28 ' 29 To obtain a proper fitting function, we first an- 
alyze separately the behaviors for large and small vortex 



separations. For large separation, using the asymptotic 
form of the modified Bessel functions, Eq. (ITU1) can be 
rewritten as 



fl{d -> oo) = 7 tH (Se~ d - y/Jk-^) 



(15) 



where 7 and <5 are fitting parameters. For small separa- 
tion, our results show that a power function of d describes 
the force satisfactorily, i.e. 



fi(d ->•()) =pd q , 



(16) 



with p and q as fitting parameters. Two examples of 
such fittings are shown in Fig. [2] for \i = 0.6 (bottom) 
and /i = 1.7 (top). Notice that the parameter q depends 
weakly on /j,, exhibiting values between « 2.7 and ~ 2.8 
for all values of fi considered in the V-V case. 

Following the established behavior of the interaction 
in the limiting cases, we propose a single function which 
has the above limits as limiting behaviors: 
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FIG. 2: (Color online) Comparison between the V-V interac- 
tion force as a function of the separation d obtained by the 
numerical method (symbols) and by the asymptotic expres- 
sions (curves) from Eqs. (|15II16[) , for fi = 1.7 (a) and fi — 0.6 
(b). The forces are plotted on a logio scale and the values of 
the fitting parameters p, q, 7 and 8 are given in each panel. 
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FIG. 3: (Color online) (a) Comparison between the V-V in- 
teraction force as a function of the separation d obtained by 
the numerical method (symbols) and by the fitting function 
(curves) given by Eq. (|17p . The solid (dashed) curves and the 
open triangles (squares) are the results for [i — 0.6 (1.7). (b-c) 
The results for each [i are plotted separately in a log-scale, to 
emphasize the maximal deviation of the fitting function from 
the data. 



where r\i (i = 1 - 4) are four fitting parameters. 

Fig. [3] shows the fitting obtained with Eq. (|T7|) for 
the same values of \i presented in Fig. [2] The fitting 
is not ideal for d < A, where the force becomes very 
small. Nevertheless, we found that the fitting error is 
lower than 1%. Please note that the V-V interaction 
potential, which is an integrated force, will be even more 
accurate. 

The values of the four fitting parameters are given in 
Table I, for ^ from 0.2 to 2.5. Notice that the estimated 
variance v increases with /i, and thus Eq. (|17[) should 
not be used in the extreme type-II case. Nevertheless, 
in the extreme type-II case the critical separation d c ap- 
proaches zero, as mentioned earlier, and consequently, 
the short range part of the V-V interaction force will 
not be important in such a situation. Hence, the asymp- 
totic expression £l(d) = foKi (d) frequently used in the 
literature) 30 ' 31 which can be obtained from Eq. (fTU)) by 
making /1 — > 00, is expected to provide a good description 
of the V-V interaction force in extreme type-II situations. 

We next attempt to find an analytical expression for 



TABLE I: Fitting parameters r\i and the estimated variance v for Eq. (|17[) in the V-V case, for different values of /j, 
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the fitting parameters as function of fi. Their depen- 
dence on [i is shown in Fig. 01 Three different ranges 
of n, delimited by vertical dotted lines in Fig. IU can be 
distinguished. The physical reason for the existence of 
three different behaviors of the parameters r\i as a func- 
tion of \x is the following: for type-I (/x < 1), as we ex- 
plained earlier, there must be a regime where the size of 
the attractive force peak increases with \x and another 
region where it decreases with \x. This defines the ranges 
1 (fj, < 0.5) and 2 (0.5 < fj, < 1), respectively. Range 3 is 
then the type-II regime, for fj, > 1, where the interaction 
force is repulsive. The functions ?fc(/z) in Fig. @] were 
fitted as 



»7i(m) 



= Bi( M Cl +A 1 ) 



(18a) 



TABLE II: Fitting parameters in Eqs. 
case, for three different ranges of \i. 



£[8] a-d) for the V-V 
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r, 3 {n)=A 3 + B 3 e c ^ 



rn(jj) = At + Bifi 



C'l 



(18b) 



(18c) 
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with different parameters Ai , Bi and Ci for each range 
listed in Table II. These fitting functions for r)i(fi) are 
shown as solid curves in Fig. [4] Notice that the pa- 
rameter ?74 must satisfy the condition 774 < ^JJL (> y/Jt) 
in the type-I (type-II) case, otherwise the difference be- 
tween the exponential terms in Eq. (|17[) would exhibit a 



sign change for small separations, leading to a spurious 
repulsive (attractive) region in this case. In the type-II 
case, this condition leads to 774 (/x) ~ ^[Jl as the best value 
for this fitting parameter. 

It is important to point out that the results obtained 
for r/2 are not the same as the values of q in Eq. (|16[) for 
the power law at small separations, which, as mentioned 
earlier, are between ss 2.7 and ~ 2.8. This is reasonable, 
because the exponential terms in Eq. (fTT)) still play a 
role in the small d limit of this expression, thus, the pa- 
rameter r) 2 must assume a value that is different from q 
in order to compensate these terms. The values obtained 
for 774, which is the parameter controlling the large d 



range of Eq. (jT0|) . are also not the same as the values 
obtained when one uses the asymptotics of each single 
vortex to find the parameters 7} and 72 in Eq. JTOj). 
Actually, for fj, > 1, we found 774 rj «//Z, which is equiva- 
lent as making 7} = 7 2 . As fi increases, the difference 
between 7^ and 72 is amplified^, leading to a higher 
variance v for large fi, as shown in Table I. Even so, this 
choice of 774 conveniently leads to a function which decays 
exponentially for large separation d, as expected for V-V 
interactions in bulk superconductors, and which exhibits 
no sign change at small separations. Of course, the fitting 
function Eq. (fT7|) can be improved to provide a better 
fitting of the large separation part and to reproduce a 
perfect power law for small separations, but this would 
require more fitting parameters and very complicated ex- 
pressions. Equation (fTT)l is simple and still accurate for 
< fi < 2.5, as verified by the small variances v < 10~ 6 
in Table I and by the comparison with the numerical re- 
sults in Fig. [3j 
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FIG. 4: (Color online) Fitting parameters (symbols) in Eq. 
(|17|l as a function of fi — \/2k for the V-V case. The curves are 
the fitting functions given by Eqs. (1 181) for r/i(fi) (i — 1 — 4), 
in three different regions: < /z < 0.5, 0.5 < /x < 1 and 
1 < fi. The inset shows a magnification of the results for 774 
at large fi. 




FIG. 5: (Color online) Numerically obtained V-GV interac- 
tion force fi as a function of the separation d between a vortex 
and a double-quantized vortex, for several values of ft — \/2k, 
in the (a) type-II and (b) type-I regimes, (c) Critical separa- 
tion d c (open squares, right scale) and extremum fi ma:E (open 
triangles, left scale), which correspond respectively to the po- 
sition and amplitude of the peak in the force, as a function of 
the GL parameter. The fitting functions for d c and Q. m ax are 
plotted by the solid curves. 



B. Vortex-giant vortex interaction 

As observed in Fig. [TJ the interaction force between 
two vortices shows a maximum at some critical separa- 
tion d c and decays to zero for both very large and very 
small separations. The former is reasonable, since the in- 
teraction between vortices is expected to weaken as they 
are placed further from each other. The latter is due to 
the formation of a giant vortex state: when two vortices 
of winding numbers e.g. n\ = 1 and n-i = 1 are put close 
to each other, they coalesce, forming a n = n\ + n-z = 2 
giant vortex.— ~— In the absence of lateral confinement, a 
giant vortex is a stable (unstable) state in type-I (type- 
II) systems and can interact as such with other vortices, 
and this motivated us to investigate the interaction force 
between a vortex and a giant vortex. 

The V-GV interaction force is shown in Fig. [5] as a 
function of the distance between them, for several val- 
ues of fi in the type-II (a) and type-I (b) regimes. The 
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FIG. 6: (Color online) Comparison between the V-GV inter- 
action force as a function of the separation d, obtained by the 
numerical method (symbols) and by the asymptotic expres- 
sions (curves) from Eqs. (| 151 1 16f) . for fj, = 1.7 (a) and fj, — 0.6 
(b). The forces are plotted on a logw scale and the values of 
the fitting parameters p, q, 7 and 8 are given in each panel. 



behavior of the curves is quite similar to those in Fig. 
[T] for the V-V case, but with different amplitudes and 
critical separations. The critical separation d c , beyond 
which the vortices start to coalesce, obtained numeri- 
cally for the V-GV interaction is shown as a function of 
the GL parameter \x in Fig. [5](c), along with its fitting 
function d c = 25.043(1 + 6.632^T - 8862 (with estimated 
variance v k, 2%). Notice that the critical separation for 
the V-GV interaction is always larger than the one for 
the V-V case, because the giant vortex has a larger core 
in comparison to a n = 1 vortex. Nevertheless, the fit- 
ting function shows that the smallest critical separation 
for the V-GV interaction force, which would be obtained 
in an extreme type-II regime, is also d c (fi —¥ 00) = 0, 
as in the V-V case. The behavior of the extremum of 
the force peak VL max as a function of /1, shown as open 
triangles in Fig. [5jc) , is similar to that found for the V- 
V case, with the amplitude approaching zero for jj, — > 
and jj, — > 1, and increasing monotonically for /i increasing 
above 1. The extremum of the force peak can be fitted 
to fl max = 0. 1709 n(fi - 1)/(1 + 1.854/x) 06087 , which is 
shown by the solid curve in Fig. [SJc). 



In Sec. II, we analytically found that Eq. ([TO]) remains 
valid for the asymptotic V-GV interactions, simply by 
choosing n\ = 1, ni = 2 and changing the parameters 

(i) (i) 

7 j and 72 accordingly. Moreover, our results show 
that the force in the small separation limit in this case 
can still be well described by a power function of the 
separation d. The fitting of the force for the small and 
large separation limits, given by Eqs. (115U16[) . are shown 
in Fig. [5]for the V-GV interaction, agreeing well with the 
numerical results. This suggests that the fitting function 
given by Eq. p7|) can not only be used for the V-V, but 
also for the V-GV interaction force. 

The four fitting parameters r\i (i = 1 - 4) found for 
each value of /j, in the V-GV case are shown in Table III, 
for /i from 0.2 to 2.5. As in the V-V case, the estimated 
variance v increases for fj, > 1, hence the offered function 
is expected to fail in the extreme type-II case. As an 
example, the V-GV interaction force for /1 = 1.7 and 
0.6 is shown in Fig. [7] as a function of the vortex-giant 
vortex separation d, along with the fitting curves given 
by Eq. (|17[) with the corresponding parameters in Table 
III. Although the estimated variances for these cases are 
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FIG. 7: (Color online) (a) Comparison between the V-GV 
interaction force as a function of the separation d obtained by 
the numerical method (symbols) and by the fitting function 
(curves) given by Eq. (|17p . The solid (dashed) curves and the 
open triangles (squares) are the data for fi = 0.6 (1.7). The 
results for each /i are also plotted separately on a log-scale: 
(b) n= 1.7; (c) [1 = 0.6. 
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smaller than 10 -5 , it can be seen in the log-plot in Fig. 
[7](b) and (c) that for small separation the fitting function 
is less accurate as compared to the V-V case shown in Fig. 
[31 where the variances are lower than 10 -7 . Nevertheless, 
in the low d region the force is small and consequently 
the deviation in the force will also be small. 

The dependence of the fitting parameters on /i is shown 
in Fig. E] (squares) for the case of a n\ = 1 and ni = 2 
V-GV interaction, where the fitting to the data is also 
shown (curves). Once more, the three different behaviors 
of r/i as a function of /j are observed, with data fitted with 
different parameters A i: Bi and d, using Eqs. fI8k-d). 
The reason for the existence of three identifiable behav- 
iors of the parameters r\i as a function of /z is the same as 
in the V-V case. Moreover, for /j > 1 (type-II regime), 
we found 774 (/i) rj «//J, similar to the V-V case. The pa- 
rameters Ai, Bi and Ct for the V-GV case, for each range 
of fi, are given in Table IV. 
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C. The three vortex problem 

Up to now, we considered only the two-body interac- 
tion of vortices. In the study of the dynamics of many 
vortices, one generally considers the sum of pairwise in- 
teractions. In this sense, the force acting on vortex i in 
a system with many vortices forming a certain configu- 
ration is given b y 30 ' 31 



a 



5> 



i\)h. 



(19) 



is the position of the vortex i and 



«ij 



where fl 

(fi — fj) j\fi —fj\. In such a model, the V-V interaction 
force is derived from the interaction potential between a 
pair of vortices, and is usually taken as fl(\fl — fj\) = 
foKi (\fl — rJI/A), where /q is a constant. As mentioned 
in Sec. IV A, this corresponds to the extreme type-II sit- 
uation (/i> 1). For intermediate values of [i, one should 
consider both Bessel functions in Eq. PH)) . which is not 
convenient since, as mentioned in Sec. II, this expres- 
sion diverges for small V-V separation and does not take 
into account neither the vortex core deformations, nor 
the formation of giant vortices. Hence, using Eq. (|T7|) 
for the vortex pair interaction force, with the parame- 
ters given by Table I or by Eqs. (ITSa -d), would be an 
easy way to take these features into account and avoid 
the small separation divergence. Although this solves the 
problem in the standard simulations of vortex dynamics, 
which involves only pairwise interactions, here we will 
show when such an approach breaks down and take the 
interaction between three vortices placed in the vertices 
of an equilateral triangle as an example. 

Let us consider three vortices placed in the positions 
ri = (-d/2, -V3d/A), r 2 = (d/2,- v / 3d/4) and r 3 = 
(0, V3d/4), forming a triangle of side d. The three vortex 
ansatz is * = e mi6l e m2e2 e m3e;s f(x,y), where we control 
the vorticity ni, n% and 71.3 of each vortex. For the present 
example, we consider three singly quantized vortices n\ — 
n 2 = n 3 — 1- Following the procedure described in Sec. 
IIII1 we obtain the Euler-Lagrange equations for the three 
vortex problem, where the first equation is similar to Eq. 
(|13a|i . but X and Y contains three terms 
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and the second equation is 
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FIG. 8: (Color online) The fitting coefficients 

f?i=i-4 (symbols) in Eq. (I17|) as a function of n — \f2n 
for the V-GV case. The curves present the rji(fi) fittings 
given by Eqs. (|18b -d) for three different ranges, defined in 
the text as range 1 (0 < fj, < 0.5), range 2 (0.5 < /j, < 1) and 
range 3 (1 < fi). 



Similarly to the case of V-V and V-GV interactions, 
we solved the Euler-Lagrange equations for three vortices 
numerically, by means of a finite-difference scheme and a 
relaxation method. The results obtained for the interac- 
tion energy as a function of the V-V distance, or equiv- 
alent^, the side d of the triangle, are shown as dashed 
curves in Fig. [S]for two values of the GL parameter, /i = 
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TABLE III: Fitting parameters r\i and estimated variance v for Eq. i|17[l in the V-GV case (ni = 1 and na = 2), for [i between 
0.2 and 2.5. 
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1.095 

1.140 

1.183 

1.225 
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0.8 (type-I) and 1.7 (type-II). As we consider the same 
distance d between each pair of vortices forming the tri- 
angle, the standard procedure for the many- vortex prob- 
lem, which considers only pair interactions, predicts an 
interaction energy E int (d) = 3Ef^ r (d), where Ef^ r (d) 
is the interaction energy for each V-V pair. This energy 
is shown by the solid curves in Fig. [9l where good agree- 
ment with the results obtained from the three vortex 
ansatz is observed only for larger separations d, whereas 
for smaller separations the energies predicted by the pair 
interaction model are clearly overestimated. This result 
is a manifestation of the importance of the vortex de- 



TABLE IV: Fitting parameters in Eqs. (HHJi-d) for the V-GV 
case (ni = 1 and n-z — 2), for three different ranges of \x. 
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FIG. 9: (Color online) Interaction energy for three vortices 
placed in the vertices of an equilateral triangle as a function 
of its side d, obtained by the three vortices ansatz (dashed) 
and by considering only interaction between pairs (solid), i.e. 
E int (d) = 3Ef™ r {d), for /j, = 0.8 and 1.7. 



formations for small V-V separation: the pairs interac- 
tion model simply does not account for giant vortex de- 
formations with three vortices. As a consequence, this 
model overestimates the energy. This is illustrated in 
Fig. 1101 where the amplitude of the order parameter 
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for the two (a, c) and three (b, d) interacting vortices 
is shown for /.i = 0.8 (left panels) and 1.7 (right pan- 
els) at V-V separations d = 3.2A and 1.8A, respectively. 
In the case of three vortices, we observe that each vor- 
tex is deformed towards the center of the vortex cluster. 
Such a deformation, which is found as the lowest energy 
state of the three vortex system, cannot be obtained by 
a model consisting only of interactions between pairs of 
vortices. Nonetheless, in the extreme type-II cases stud- 
ied in the literatur o 30 ' 31 , the critical V-V separation d c 
where the vortices start to coalesce approaches zero (as 
demonstrated in the previous subsection) and the agree- 
ment between the results obtained by the pairwise model 
and by the three vortex ansatz is expected to improve at 
even smaller separations. 

In the many- vortex problem, the vortices can approach 
each other in a combinatorially large number of ways, 
and the study of three vortices in a triangular geometry 
presented in this subsection is a very specific case. Nev- 
ertheless, this example illustrates in a simple way that, 
apart from the extreme type-II case, an exact description 
of the many- vortex dynamics is a very difficult task. The 
pair potential, even when taking V-V deformations into 
account, still provides only an approximate description 
of the problem, as the deformations due to the presence 
of all the other vortices are not included in the model. 
In this sense, the expressions proposed in the present 
work for the V-V and V-GV interaction forces provide an 
important improvement on the well-known expressions 
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FIG. 10: (Color online) Amplitude of the order parameter for 
the cases indicated in Fig. |9j namely, for /j, — 0.8 (1.7) and 
d = 3.2A (1.8A), obtained for the two (a, c) and three (b, d) 
interacting vortices. 



to(\ri-fj\) = foKi {\n ~ fj\/\) and Eq. dTUJ) , by includ- 
ing the deformations and the merger of vortices, but a 
molecular dynamics study of many vortices using these 
expressions is still not an ideal description of a system 
with comparable length scales £ and A. 



D. Vortex-antivortex interaction 

In the previous subsections, we showed that when two 
vortices or a vortex and a giant vortex are brought close 
to each other, they merge forming a single giant vortex 
state with vorticity n = ni + n®, and in the limit of small 
separation the V-V or V-GV forces are very weak. Con- 
versely, as discussed in Sec. II, a vortex and an antivortex 
attract and annihilate, both in type-I and type-II super- 
conductors. In what follows, the behavior of the force for 
the vortex-antivortex (V-AV) interaction as a function of 
the V-AV separation is studied in greater detail. 

Indeed, the V-AV interaction is quite different from the 
interaction observed in the V-V and V-GV cases studied 
in previous subsections. The V-AV interaction energy (a) 
and force (b) are shown in Fig. [TTJ for two values of the 
GL parameter, /j, = 0.6 (triangles) and 1.7 (squares). As 
discussed previously in Sec. II, the V-AV interaction is 
always attractive, for any value of \x. However, at some 
critical V-AV separation de, the solution with well de- 
fined super-currents around each vortex and antivortex 
ceases to be the lowest energy state of the system. A so- 
lution with lower energy exhibits a strong suppression of 
the amplitude of the order parameter and super-current 
in the region between the vortex and the antivortex, and 
represents the ground state for small V-AV distances. A 
hysteresis is observed in the vicinity of the critical sep- 
aration dE, as shown in Fig. UTT a). These results re- 
semble those obtained by Priour and Fertig 3 ^ in the case 
of a vortex placed close to an artificial defect. A sup- 
pressed amplitude of the order parameter is also observed 
by Sardella et at 33 in the dynamics of V-AV annihila- 
tion in a square mesoscopic superconducting cylinder, for 
small V-AV separation. The absolute value of the force 
is shown in Fig. Illf b) on logio scale, where two different 
behaviors, for separations d smaller and larger than ds, 
are clearly observed. 

The dependence of the numerically obtained critical 
separation d^ for the V-AV interaction on the GL pa- 
rameter [i is illustrated as the squares in Fig. Illf c). 
and can be fitted to a function similar to the one used 
for the critical separations in the V-V and V-GV cases, 
given by d E = 0.337 + 31.249(1 + 10.264/i)-° 6855 (with 
estimated variance v « 0.4%), which is shown as a 
solid curve in Fig. 10(c). The difference is that the 
fi — > oo limit is now finite, while previously it was zero. 
Although the solution with well defined super-currents 
around each vortex and antivortex is not the lowest en- 
ergy state for d < d^, it is still a stable state in the 
vicinity of this point, and becomes unstable only for 
d < dA- The dependence of dA on the GL parameter 
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FIG. 11: (Color online) V-AV interaction (a) energy and (b) 
force (absolute value) as a function of the separation d, for 
H = 0.6 (triangles) and 1.7 (squares). The full (open) symbols 
are the results obtained in the numerical relaxation procedure 
by gradually increasing (decreasing) d from to 15 A (from 15 
A to 0). A hysteresis is observed around a critical separation 
d_E, as indicated by the arrows, and the solution represented 
by open symbols is stable only for d > cLa- (c) Numerically 
obtained critical separations d_e (squares) and (La (triangles) 
as a function of the GL parameter n, along with their fitting 
functions (curve). 



fi is shown by the triangles in Fig. [TTTc) and can be 
fitted by d A = 0.337 + 12.222(1 + 2.461/i) 07931 (with 
estimated variance v k, 0.6%). 

Figure Q2] shows the distribution of the super-current 
J = V x V x A and the amplitude of the order parameter 
along the direction of the V-AV approach (y = axis) 
for different values of the V-AV separation. The criti- 
cal separations for /i = 0.6 and 1.7 are d,E = 8.6A and 
4.5A, respectively, and the values of V-AV separation in 
Fig. [T5] are chosen as d = 9.2A > dE for /i = 0.6 and 
5. 2 A > dE for [i = 1.7. Notice that for each of these 
separations, we can find two solutions with different en- 
ergies. The black (gray) curves in Fig. [T^are related to 
the open (full) symbols in Fig. [TlT b). When the V-AV 
separation is large, the currents around the vortex and 



the antivortex present well defined peaks at some dis- 
tance which depends on /i. As the vortex and antivortex 
are placed closer, their super-currents superimpose in the 
region between them, as observed in the black curves in 
Fig. [12] (top). The black curves in Fig. [12] (bottom) 
show that the amplitude of the order parameter in these 
solutions has zeros at each vortex and antivortex posi- 
tion and reaches sa 1 in the region between them. In the 
vicinity of dE, for V-AV separation d > dE, there is a 
higher energy state (see full gray symbols in Fig[TTTa)) 
with strongly suppressed super-current and amplitude of 
the order parameter in the region between the vortex and 
antivortex, which is shown by the gray curves in Fig. 1121 
For d < ds, the solution represented by black lines in 
Fig. [T2] is no longer the lowest energy state, as shown in 
Fig- rTTT a). and becomes unstable as the V-AV separation 
is reduced at d < d A , while the solution with suppressed 
current and order parameter between the vortex and an- 
tivortex, shown by the gray curves, becomes the lowest 
energy state for d < ds and the only stable solution for 
d < d A - 

The suppressed order parameter in the region between 
vortices observed in the only stable solution for d < d A 
suggests that a vortex and an antivortex cannot coexist 
at these distances, unless somehow pinned, in which case 
this string solution is formed. This is reasonable, since at 



n = 0.6 



H=1.7 
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FIG. 12: Super-current (top) and amplitude of the order pa- 
rameter (bottom) along the direction of the V-AV approach 
x, for the V-AV separations indicated by arrows in Fig. 10 
(b), namely, d = 9.2A (5.2 A), for /i = 0.6 (1.7). Black (gray) 
curves refer to the states represented by open (full) symbols 
in Fig. [11] (b). 
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TABLE V: Fitting parameters Aj and estimated variance v 
for Eq. (f2"Tj) in the V-AV case (m = 1 and n 2 = —1), for 
d > ds- 



/' 



Ai 



i/(xlO" 



0.3 


156.948 


0.4203 


0.918 


0.4 


67.315 


0.7419 


0.10 


0.5 


31.064 


1.173 


0.098 


0.6 


19.070 


1.719 


0.18 


0.7 


14.401 


2.340 


0.342 


0.8 


11.990 


2.925 


0.56 


0.9 


6.499 


5.060 


3.06 


1.0 


6.357 


6.357 


5.01 


1.1 


5.320 


9.170 


5.0 


1.2 


4.159 


14.269 


1.21 


1.3 


4.254 


16.759 


1.21 


1.4 


4.125 


20.640 


1.10 


1.5 


4.039 


23.58 


0.75 


1.6 


3.775 


33.662 


1.20 


1.7 


3.632 


43.173 


1.40 


1.8 


3.542 


51.251 


1.03 


1.9 


3.422 


66.755 


1.23 


2.0 


3.315 


87.448 


1.38 


2.1 


3.226 


113.229 


1.54 


2.2 


3.162 


138.472 


1.25 


2.3 


3.101 


170.634 


1.08 


2.4 


3.037 


220.086 


1.13 


2.5 


2.983 


277.742 


1.08 



short distances the fields of the vortex and the antivortex 
compensate each other, and the flux quantization as an 
essential property of a(n) (anti)vortex is lost. Notice that 
this is different from the case of two merging vortices, 
which can coexist at short distances, deform and interact 
as described in previous sections, since the flux quantiza- 
tion of the V-V pair is preserved even at small V-V sepa- 
rations. The string formation goes beyond simulations of 
V-AV dynamics, since in this case the vortex-antivortex 
pair is no longer well defined by their surrounding super- 
current and order parameter. For molecular dynamics 
studies of the V-AV motion, one should consider the crit- 
ical separation cIa as the separation where the V-AV pair 
annihilates (see, e.g. Refi^i). 

Due to the peculiar behavior found for the V-AV force 
as a function of the separation d, which is discontinuous 
at d,E, it is not possible to find a single fitting function 
describing the force for both the d > d,B and d < dE 
regimes, as we were able to do for the V-V and V-GV 
forces. On the other hand, as discussed in the previous 
section, the V-AV interaction force at large distances d 
can be described by a combination of Bessel functions, 
given by Eq. (jlOp . which can be rewritten as 



Vt{d) = -Ai-K^d) - A 2 Ki(nd), 



(21) 



where Ai and A2 are fitting parameters. We have fit- 
ted our numerically obtained V-AV interaction force for 
d > dE using Eq. (j2"Tj) . and a list of fitting parameters 
for the GL parameter /x ranging from 0.3 to 2.5 is given 
in Table V. A list of such fitting parameters can also 
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FIG. 13: (Color online) Fitting functions (curves) for the pa- 
rameters Ai and A2 in Eq. (|21[) . for the V-AV interaction 
force at d > de, as a function of the GL parameter n = \/2k. 
The data from Table [V] are shown as symbols for comparison. 



be found in Ref.— , where the relation between our fit- 
ting parameters and the parameters q and m of the cited 
work is Ai = m 2 /2tt 2 and A 2 = [iq 2 /2tt 2 . Following the 
same procedure of previous subsections, we propose fit- 
ting functions for these parameters as a function of the 
GL parameter /x, 



Ai(/x) = 2.879 + 3.415ii 



-3.166 



and 



A 2 (aO = M-0.2258+ 1.044e 1866 ^), 



(22a) 



(22b) 



which are plotted as solid lines in Fig. Q2] along with the 
data of Table [V] (symbols). 

Notice that for the V-AV interaction, we did not find 
different behaviors in different ranges of /i, as observed 
for the V-V and V-GV cases, since the V-AV interac- 
tion is always attractive and becomes only stronger as 
fi increases from zero, instead of exhibiting zero force at 
/1 = 1 and becoming repulsive for fi > 1, as observed 
for V-V and V-GV interactions. Therefore, substituting 
Eqs. (|22|) into Eq. (|2"Tj) yields a single expression 



Q(d) = -(2.879 + 3.415// 



-3.166 



)KM 



+ (0.2258 - 1.044e 1 - 866 '>lfi(/id), 



(23) 



which is expected to provide an accurate description of 
the V-AV interaction force, at separations d > dE, for 
any value of \x. 



V. CONCLUSIONS 

We presented a theoretical study of the interaction 
between vortices in bulk superconductors within the 
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Ginzburg-Landau theory. An analytical study of the 
asymptotic behavior of the vortex-vortex interaction 
shows that a combination of first order modified Bessel 
functions of the second kind describes the behavior of the 
numerically obtained forces for large vortex-vortex sep- 
aration. At small distances, the fitting curves for V-V 
and V-GV interactions show that the force in this region 
behaves as a power function of the separation between 
vortices. We proposed a fitting function that combines 
both limiting behaviors, namely, the power law for small 
distances and the modified Bessel function behavior for 
large distances. This function, given by Eq. (fTT|). gives 
fairly accurate fitting of the interaction force for any value 
of /i, even in the type-I regime. It depends on four fit- 
ting parameters, which can be obtained for any value of 
/i either by interpolating our data presented in Table I 
(Table III), for vortex- vortex (vortex-giant vortex) inter- 
actions, or by using Eqs. (fT51 a-d) with the parameters 
presented in Table II (Table IV). 

Our analytical study of the V-AV interaction shows 
that the V-AV interaction force is attractive for any value 
of the GL parameter //, which is confirmed by our nu- 
merical results and contradicts the conjecture proposed 
in previous worksi£~— which implies that the V-AV in- 
teraction force is repulsive for \x < 1 (type-I). For large 
V-AV separation d, the interaction force decays with d as 
a combination of modified Bessel functions. However, for 
d smaller than a critical separation dE, the conventional 
V-AV pair is no longer the lowest energy state. Instead, 
the lowest energy state exhibits a strong suppression of 
the super-current and amplitude of the order parameter 
in the region between the vortex and antivortex, which 
results in a different behavior of the force as a function of 
d in this case and, as a consequence, the V-AV interaction 
force is discontinuous at dE- Furthermore, the conven- 



tional V-AV pair becomes unstable for d lower than the 
separation d^ which is interpreted as the V-AV annihi- 
lation point. We fitted the interaction force for V-AV 
separations d > dE by Eq. ([2"Tj) and proposed an ap- 
proximate analytical expression for the V-AV interaction 
force at these separations, given by Eq. (f23"l) . which is 
valid for any value of [x. 

The fitting functions for the V-V, V-GV and V-AV 
force given in this work will be useful, for instance, for the 
study of bulk and mesoscopic systems consisting of many 
vortices using molecular dynamics techniques. We nev- 
ertheless remark that, although deformations are taken 
into account in the interaction force between two vortices 
in this work, the deformations in a many vortices system 
are expected to be more complex. Hence, the molecular 
dynamics study of many vortices, even with the improved 
expressions for the interaction force provided in this pa- 
per, is still an approximate description of the system. As 
a method, the derivation and handling of the differen- 
tial equations describing the interaction between vortices 
presented in this work can be further adapted to describe 
such interactions in e. g. two-band superconductors, or 
hybrid systems comprising different superconducting ma- 
terials. 
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